#####Load packages

#load packages

#library(tidyverse) #there is an error with tidyverse
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(devtools)
## Warning: package 'devtools' was built under R version 4.1.3
## Loading required package: usethis
## Warning: package 'usethis' was built under R version 4.1.3
library(stats)
library(ordinal)
## Warning: package 'ordinal' was built under R version 4.1.3
## 
## Attaching package: 'ordinal'
## The following object is masked from 'package:dplyr':
## 
##     slice
library(MASS)
## 
## Attaching package: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
library(ggeffects)
## Warning: package 'ggeffects' was built under R version 4.1.3
library(effects)
## Warning: package 'effects' was built under R version 4.1.3
## Loading required package: carData
## Warning: package 'carData' was built under R version 4.1.3
## lattice theme set by effectsTheme()
## See ?effectsTheme for details.
library(scales)
## Warning: package 'scales' was built under R version 4.1.3
library(car) #for regression diagnostics
## Warning: package 'car' was built under R version 4.1.3
## 
## Attaching package: 'car'
## The following object is masked from 'package:dplyr':
## 
##     recode
library(broom) #for tidy output
## Warning: package 'broom' was built under R version 4.1.3
library(ggfortify) #for model diagnostics
## Warning: package 'ggfortify' was built under R version 4.1.3
#library(sjPlot) #for outputs
library(knitr) #for kable
## Warning: package 'knitr' was built under R version 4.1.3
library(effects) #for partial effects plots
library(MuMIn) #for interfacing with STAN
library(emmeans) #for estimating marginal means
## Warning: package 'emmeans' was built under R version 4.1.3
## 
## Attaching package: 'emmeans'
## The following object is masked from 'package:devtools':
## 
##     test
library(ggeffects) #for partial effects plots
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(DHARMa) #for residual diagnostics plots
## Warning: package 'DHARMa' was built under R version 4.1.3
## This is DHARMa 0.4.6. For overview type '?DHARMa'. For recent changes, type news(package = 'DHARMa')
library(modelr) #for auxillary modelling functions
## Warning: package 'modelr' was built under R version 4.1.3
## 
## Attaching package: 'modelr'
## The following object is masked from 'package:broom':
## 
##     bootstrap
## The following object is masked from 'package:ggeffects':
## 
##     data_grid
library(performance) #for residuals diagnostics
## Warning: package 'performance' was built under R version 4.1.3
## 
## Attaching package: 'performance'
## The following objects are masked from 'package:modelr':
## 
##     mae, mse, rmse
library(see) #for plotting residuals
## Warning: package 'see' was built under R version 4.1.3
library(patchwork) #for grids of plots
## Warning: package 'patchwork' was built under R version 4.1.3
## 
## Attaching package: 'patchwork'
## The following object is masked from 'package:MASS':
## 
##     area
library(brms)
## Warning: package 'brms' was built under R version 4.1.3
## Loading required package: Rcpp
## Warning: package 'Rcpp' was built under R version 4.1.3
## Loading 'brms' package (version 2.18.0). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').
## 
## Attaching package: 'brms'
## The following object is masked from 'package:MuMIn':
## 
##     loo
## The following object is masked from 'package:stats':
## 
##     ar
library(bayesplot)
## Warning: package 'bayesplot' was built under R version 4.1.3
## This is bayesplot version 1.10.0
## - Online documentation and vignettes at mc-stan.org/bayesplot
## - bayesplot theme set to bayesplot::theme_default()
##    * Does _not_ affect other ggplot2 plots
##    * See ?bayesplot_theme_set for details on theme setting
library(coda)
## Warning: package 'coda' was built under R version 4.1.3
library(ggmcmc)
## Warning: package 'ggmcmc' was built under R version 4.1.3
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(broom) #for tidying outputs
library(tidybayes) #for more tidying outputs
## 
## Attaching package: 'tidybayes'
## The following objects are masked from 'package:brms':
## 
##     dstudent_t, pstudent_t, qstudent_t, rstudent_t
library(broom.mixed)#for summarising models
## Warning: package 'broom.mixed' was built under R version 4.1.3
library(knitr)

#Arranging data set

Data 2016, treatment and control communities.

Change name of the data set

data_fiji<- read.csv ("Fiji_data_notabu_2016.4.csv",
                     header=TRUE,
                    dec=".")

Transform variables into factors (District, control/treatment)

data_fiji$village<-as.factor(data_fiji$village) #village as factor

dim(data_fiji) # n=193
## [1] 193 131

COVARIATES:

AGE

#Age (nominal)

str(data_fiji$Age)
##  chr [1:193] "33" "35" "63" "69" "50" "40" "32" "34" "77" "73" "46" "60" ...
ggplot(data_fiji,aes(Age))+geom_bar()

data_fiji$Age_num<-as.numeric(as.character(data_fiji$Age)) 
## Warning: NAs introduced by coercion
hist(data_fiji$Age_num) 

MIGRANT

#Migrant status (Categorical)

str(data_fiji$Migrant)
##  int [1:193] 3 1 2 1 1 1 1 3 1 1 ...
#transform in factor
data_fiji$Migrant<-as.factor(data_fiji$Migrant)
str(data_fiji$Migrant)
##  Factor w/ 4 levels "1","2","3","4": 3 1 2 1 1 1 1 3 1 1 ...
#Combine Migrant levels 2,3 and 4

data_fiji$Migrant<-as.numeric(as.character( data_fiji$Migrant))
str(data_fiji$Migrant)
##  num [1:193] 3 1 2 1 1 1 1 3 1 1 ...
data_fiji$migrant2 <- as.factor(ifelse(data_fiji$Migrant > 1,2,1)) 
data_fiji$Gender<-as.factor(data_fiji$Gender)
#Marital status

str(data_fiji$Marital_status)
##  chr [1:193] "m" "m" "m" "s" "m" "s" "s" "s" "m" "m" "m" "m" "m" "s" "w" ...
data_fiji$Marital_status<-as.factor(data_fiji$Marital_status)
#Education 

str(data_fiji$Education)
##  chr [1:193] "Tertiary" "Secondary " "Secondary" "Primary" "Secondary" ...
data_fiji$Education<-as.factor(data_fiji$Education)
str(data_fiji$Education)
##  Factor w/ 25 levels "0","Class 2",..: 25 21 18 16 18 18 18 18 16 16 ...
data_fiji$Education_num<-as.factor(ifelse(data_fiji$Education == "Primary",1,
                                                   ifelse(data_fiji$Education=="Secondary",2,
                                                          ifelse(data_fiji$Education=="Tertiary",3,0))))


data_fiji$Education_num<-as.numeric(as.character(data_fiji$Education_num))

str(data_fiji$Education_num)
##  num [1:193] 3 0 2 1 2 2 2 2 1 1 ...
#Transform in Primary, secondary and tertiary

data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary-F5" = "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"      "0"              "Class 2"        "class 5"       
##  [5] "Class 6"        "class 7"        "Class 7"        "class 8"       
##  [9] "Class 8"        "Form 3"         "Form 4"         "Form 5"        
## [13] "Form 6"         "Form 7"         "Intermediate"   "primary"       
## [17] "Primary"        "secondary"      "Secondary-F6"   "Secondary "    
## [21] "Secondary - F1" "Secondary - F4" "tertiary"       "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary-F6" = "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"      "0"              "Class 2"        "class 5"       
##  [5] "Class 6"        "class 7"        "Class 7"        "class 8"       
##  [9] "Class 8"        "Form 3"         "Form 4"         "Form 5"        
## [13] "Form 6"         "Form 7"         "Intermediate"   "primary"       
## [17] "Primary"        "secondary"      "Secondary "     "Secondary - F1"
## [21] "Secondary - F4" "tertiary"       "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary - F1" = "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"      "0"              "Class 2"        "class 5"       
##  [5] "Class 6"        "class 7"        "Class 7"        "class 8"       
##  [9] "Class 8"        "Form 3"         "Form 4"         "Form 5"        
## [13] "Form 6"         "Form 7"         "Intermediate"   "primary"       
## [17] "Primary"        "secondary"      "Secondary "     "Secondary - F4"
## [21] "tertiary"       "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "Secondary - F4" = "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "0"            "Class 2"      "class 5"      "Class 6"     
##  [6] "class 7"      "Class 7"      "class 8"      "Class 8"      "Form 3"      
## [11] "Form 4"       "Form 5"       "Form 6"       "Form 7"       "Intermediate"
## [16] "primary"      "Primary"      "secondary"    "Secondary "   "tertiary"    
## [21] "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "primary" = "Primary")
levels(data_fiji$Education)
##  [1] "Primary"      "Secondary"    "0"            "Class 2"      "class 5"     
##  [6] "Class 6"      "class 7"      "Class 7"      "class 8"      "Class 8"     
## [11] "Form 3"       "Form 4"       "Form 5"       "Form 6"       "Form 7"      
## [16] "Intermediate" "secondary"    "Secondary "   "tertiary"     "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "secondary"= "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "Primary"      "0"            "Class 2"      "class 5"     
##  [6] "Class 6"      "class 7"      "Class 7"      "class 8"      "Class 8"     
## [11] "Form 3"       "Form 4"       "Form 5"       "Form 6"       "Form 7"      
## [16] "Intermediate" "Secondary "   "tertiary"     "Tertiary"
data_fiji$Education<-recode_factor(data_fiji$Education, "tertiary"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Form 3"       "Form 4"       "Form 5"       "Form 6"      
## [16] "Form 7"       "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 4"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Form 3"       "Form 5"       "Form 6"       "Form 7"      
## [16] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 3"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Form 5"       "Form 6"       "Form 7"       "Intermediate"
## [16] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 5"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Form 6"       "Form 7"       "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 6"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Form 7"       "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Form 7"= "Tertiary")
levels(data_fiji$Education)
##  [1] "Tertiary"     "Secondary"    "Primary"      "0"            "Class 2"     
##  [6] "class 5"      "Class 6"      "class 7"      "Class 7"      "class 8"     
## [11] "Class 8"      "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 2"= "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 5"     
##  [6] "Class 6"      "class 7"      "Class 7"      "class 8"      "Class 8"     
## [11] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 5"= "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 5"     
##  [6] "Class 6"      "class 7"      "Class 7"      "class 8"      "Class 8"     
## [11] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 6"= "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 5"     
##  [6] "class 7"      "Class 7"      "class 8"      "Class 8"      "Intermediate"
## [11] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 7"= "Secondary")
levels(data_fiji$Education)
##  [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 5"     
##  [6] "class 7"      "class 8"      "Class 8"      "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 7"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 5"     
## [6] "class 8"      "Class 8"      "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 5"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary"    "Tertiary"     "Primary"      "0"            "class 8"     
## [6] "Class 8"      "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "class 8"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary"    "Tertiary"     "Primary"      "0"            "Class 8"     
## [6] "Intermediate" "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Class 8"= "Secondary")
levels(data_fiji$Education)
## [1] "Secondary"    "Tertiary"     "Primary"      "0"            "Intermediate"
## [6] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "0"= "None")
levels(data_fiji$Education)
## [1] "None"         "Secondary"    "Tertiary"     "Primary"      "Intermediate"
## [6] "Secondary "
data_fiji$Education<-recode_factor(data_fiji$Education, "Intermediate"= "Primary") 
levels(data_fiji$Education)
## [1] "Primary"    "None"       "Secondary"  "Tertiary"   "Secondary "
ggplot(data_fiji,aes(Education))+geom_bar()

We will combine none and primary because there are very few people with no education

data_fiji$Education_cat3<-as.factor(ifelse(data_fiji$Education == "Tertiary",3,
                                                   ifelse(data_fiji$Education=="Secondary",2,1)))

Wealth

I have calculated wealth using all assets and chose the one with loading factors higher than 0.4 (Gurney and Darling 2017. A Global Social-Ecological Systems Monitoring Framework for Coastal Fisheries Management: A Practical Monitoring Handbook)

Loading factor is how much each variable is correlated to axis 1.

calculate wealth (Material lifestyle)

#dvd
str(data_fiji$dvd)
##  chr [1:193] "y" "y" "y" "n" "n" "y" "y" "n" "y" "y" "y" "n" "n" "n" "n" ...
data_fiji$dvd<-as.factor(data_fiji$dvd)
which(is.na(data_fiji$dvd))
## integer(0)
#washing
str(data_fiji$washing)
##  chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$washing<-as.factor(data_fiji$washing)
which(is.na(data_fiji$washing))
## integer(0)
#computer
str(data_fiji$computer)
##  chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$computer<-as.factor(data_fiji$computer)
which(is.na(data_fiji$computer))
## integer(0)
#fridge
str(data_fiji$fridge)
##  chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "y" "n" "n" "n" "n" "n" "n" ...
data_fiji$fridge<-as.factor(data_fiji$fridge)
which(is.na(data_fiji$fridge))
## integer(0)
#tv
str(data_fiji$tv)
##  chr [1:193] "y" "y" "y" "n" "n" "y" "y" "y" "y" "n" "n" "n" "n" "n" "y" ...
data_fiji$tv<-as.factor(data_fiji$tv)
which(is.na(data_fiji$tv))
## integer(0)
#satellite
str(data_fiji$satellite)
##  chr [1:193] "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" "n" ...
data_fiji$satellite<-as.factor(data_fiji$satellite)
data_fiji<-data_fiji%>% 
  mutate(satellite=na_if(satellite, "no answer"))
which(is.na(data_fiji$satellite))#187 is missing
## [1] 187
#recode rest of the variables

data_fiji$dvd<-recode_factor(data_fiji$dvd,"y" ="1")
data_fiji$dvd<-recode_factor(data_fiji$dvd,"n" ="0")

data_fiji$smartphone<-recode_factor(data_fiji$smartphone,"y" ="1")
data_fiji$smartphone<-recode_factor(data_fiji$smartphone,"n" ="0")

data_fiji$washing<-recode_factor(data_fiji$washing,"y" ="1")
data_fiji$washing<-recode_factor(data_fiji$washing,"n" ="0")

data_fiji$fridge<-recode_factor(data_fiji$fridge,"y" ="1")
data_fiji$fridge<-recode_factor(data_fiji$fridge,"n" ="0")

data_fiji$tv<-recode_factor(data_fiji$tv,"y" ="1")
data_fiji$tv<-recode_factor(data_fiji$tv,"n" ="0")

data_fiji$satellite<-recode_factor(data_fiji$satellite,"y" ="1")
data_fiji$satellite<-recode_factor(data_fiji$satellite,"n" ="0")

#transform variable into numeric vars
data_fiji$dvd<-as.numeric(as.character(data_fiji$dvd))
which(is.na(data_fiji$dvd))
## integer(0)
data_fiji$washing<-as.numeric(as.character(data_fiji$washing))
which(is.na(data_fiji$washing))
## integer(0)
data_fiji$fridge<-as.numeric(as.character(data_fiji$fridge))
which(is.na(data_fiji$fridge))
## integer(0)
data_fiji$tv<-as.numeric(as.character(data_fiji$tv))
which(is.na(data_fiji$tv))
## integer(0)
data_fiji$satellite<-as.numeric(as.character(data_fiji$satellite))
which(is.na(data_fiji$satellite)) #missing 187
## [1] 187
#PCA 18 only have one missing value.

data_fiji<-data_fiji[-c(187),]
wealth_rows_pca18.2= data_fiji[, c("dvd", "washing","fridge","tv", "satellite")]
xord18.2<-prcomp(wealth_rows_pca18.2, center= TRUE, scale=TRUE)

summary(xord18.2)
## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5
## Standard deviation     1.6008 1.0004 0.7855 0.67895 0.59890
## Proportion of Variance 0.5125 0.2002 0.1234 0.09219 0.07174
## Cumulative Proportion  0.5125 0.7127 0.8361 0.92826 1.00000
biplot(xord18.2)

loadings18.2<-xord18.2$rotation #rotation is the matrix of variable loadings (columns are eigenvectors) in prcomp() function
loadings18.2 
##                  PC1        PC2         PC3         PC4         PC5
## dvd       -0.4456042  0.5510784 -0.10704930  0.23661393 -0.65597547
## washing   -0.4359240 -0.4797544  0.32211120  0.68440724  0.08739014
## fridge    -0.4779007 -0.1855376  0.52473662 -0.66155942 -0.15549067
## tv        -0.4702871  0.4877233 -0.06487272 -0.02063843  0.73233979
## satellite -0.4023035 -0.4402841 -0.77796391 -0.19368470 -0.03949986
#extract pc1


myscores.2<-xord18.2$x # "x" are the coordinates of the individuals (observations) on the principal components in prcomp()

data_fiji$wealth_msl.2<-myscores.2[,"PC1"] #add this to fiji_data and add na to the empty cells

data_fiji$wealth_msl.2
##   [1] -0.7388053 -0.7388053 -0.7388053  1.2056164  1.2056164 -0.7388053
##   [7] -0.7388053  0.2339693 -1.8918196  0.2328417  0.2328417  1.2056164
##  [13]  1.2056164  1.2056164  0.2339693  1.2056164 -0.7388053 -0.4520388
##  [19]  1.2056164  1.2056164  1.2056164 -0.7388053  1.2056164  0.2339693
##  [25]  1.2056164 -0.7388053  1.2056164  0.0526021  1.2056164  0.2328417
##  [31]  1.2056164  0.2339693  1.2056164  1.2056164  1.2056164  1.2056164
##  [37] -1.8918196  1.2056164 -0.7388053  1.2056164  1.2056164  1.2056164
##  [43]  1.2056164  1.2056164  1.2056164  1.2056164 -0.7388053  1.2056164
##  [49]  1.2056164  1.2056164 -0.9049950 -0.7388053 -0.7388053  1.2056164
##  [55]  1.2056164  1.2056164 -3.0307839  1.2056164  1.2056164  0.0526021
##  [61]  1.2056164  1.2056164  0.2328417  1.2056164  1.2056164  1.2056164
##  [67] -1.8918196  0.0666521  1.2056164  0.2339693  1.2056164  0.0526021
##  [73]  1.2056164 -1.8918196  0.0666521  1.2056164  1.2056164  1.2056164
##  [79]  0.0666521  0.0666521 -1.8777696 -0.7388053  1.2056164  0.0666521
##  [85]  1.2056164  1.2056164  1.2056164  1.2056164  1.2056164  0.2328417
##  [91]  1.2056164  0.2339693 -2.0580092  1.2056164  1.2056164  1.2056164
##  [97]  0.2339693  1.2056164  1.2056164  0.2339693 -1.8918196 -2.0580092
## [103]  1.2056164  1.2056164 -3.7156644  1.2056164  0.2339693 -2.0580092
## [109]  1.2056164 -4.6884391  1.2056164  1.2056164 -0.7388053 -0.9190450
## [115]  1.2056164 -1.0863622 -3.0307839  0.0526021 -0.9190450  0.2328417
## [121]  0.0526021 -0.7388053  1.2056164 -0.7388053 -1.8777696  1.2056164
## [127]  1.2056164  1.2056164 -1.8918196  1.2056164  0.2339693 -0.7388053
## [133]  1.2056164  1.2056164  1.2056164  1.2056164  1.2056164  1.2056164
## [139]  1.2056164  1.2056164  1.2056164 -1.8918196  1.2056164  0.2328417
## [145] -0.7388053  0.2339693 -4.6884391  0.2339693 -3.7156644 -1.0863622
## [151] -1.8777696  0.0666521 -4.6884391 -3.0307839  1.2056164 -0.7388053
## [157]  1.2056164  1.2056164 -3.0307839 -2.3964605 -3.0307839 -4.6884391
## [163] -4.6884391 -2.0591368 -2.0591368 -4.6884391 -0.7388053 -3.0307839
## [169]  0.2339693 -1.8918196  1.2056164  1.2056164  1.2056164  1.2056164
## [175] -1.8918196 -1.8918196  0.0526021  1.2056164  1.2056164  1.2056164
## [181]  0.0666521  0.2339693  1.2056164  1.2056164 -4.6884391 -0.9190450
## [187]  1.2056164  1.2056164  1.2056164  1.2056164 -4.6884391  1.2056164
str(data_fiji$wealth_msl.2)
##  num [1:192] -0.739 -0.739 -0.739 1.206 1.206 ...
data_fiji$wealth_msl.2_sc.1<-scale( data_fiji$wealth_msl.2, center = TRUE,scale=TRUE)


data_fiji$wealth_msl.2.resc<- rescale(data_fiji$wealth_msl.2, to = c(0,1)) #rescale
data_fiji$wealth2<-data_fiji$wealth_msl.2.resc
#gender (Gender) 
str(data_fiji$Gender)
##  Factor w/ 2 levels "f","m": 1 2 2 2 2 2 1 1 2 2 ...
#marital status
str(data_fiji$Marital_status)
##  Factor w/ 3 levels "m","s","w": 1 1 1 2 1 2 2 2 1 1 ...
#age (Age)
str(data_fiji$Age_num)
##  num [1:192] 33 35 63 69 50 40 32 34 77 73 ...
#migrant status (migrant2)
str(data_fiji$migrant2)
##  Factor w/ 2 levels "1","2": 2 1 2 1 1 1 1 2 1 1 ...
#education (Education)
str(data_fiji$Education)
##  Factor w/ 5 levels "Primary","None",..: 4 5 3 1 3 3 3 3 1 1 ...
str(data_fiji$Education_cat3)
##  Factor w/ 3 levels "1","2","3": 3 1 2 1 2 2 2 2 1 1 ...
#wealth
str(data_fiji$wealth_msl.2)
##  num [1:192] -0.739 -0.739 -0.739 1.206 1.206 ...
str(data_fiji$wealth2)
##  num [1:192] 0.67 0.67 0.67 1 1 ...

Standardize continuous variables

# Scale all "numerical" variables from 0 to 1
#migrant: binomial


data_fiji$wealth2_sd <- scale(data_fiji$wealth2)*0.5
data_fiji$Age_num_sd <- scale(data_fiji$Age_num)*0.5

RESPONSE VARIABLES

I am going to create a data set for procedural equity and another one for distributional equity to avoid the elimnations of nas in procedural equity resulting from nas occuring only in disributional equity and viceversa.

DISTRIBUTIONAL EQUITY

ggplot(data_fiji, aes(Distr_Eq_Impact))+geom_bar()

Eliminate NA and don’t knows (including 6), and noanswer

data_fiji_dist <- data_fiji[data_fiji$Distr_Eq_Impact!=6,]

data_fiji_dist <- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="dontknow",]

data_fiji_dist<- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="noanswer",]

which(is.na(data_fiji_dist$Distr_Eq_Impact))
## integer(0)
data_fiji_dist <- data_fiji_dist[data_fiji_dist$Distr_Eq_Impact!="na",]
ggplot(data_fiji_dist, aes(Distr_Eq_Impact))+geom_bar()

ggplot(data_fiji_dist, aes(Distr_Eq_Impact,fill=Gender))+geom_bar()+theme_classic()

table<-xtabs(~Distr_Eq_Impact+migrant2+Gender, data=data_fiji_dist)
ftable(table)
##                          Gender  f  m
## Distr_Eq_Impact migrant2             
## 1               1                2  2
##                 2                1  1
## 2               1                6  5
##                 2                7  6
## 3               1                8  2
##                 2                3  1
## 4               1               19 36
##                 2               17  4
## 5               1               12 19
##                 2                5  4
table2<-ftable(table)
dim(data_fiji_dist)#number of observations
## [1] 160 141
which(is.na(data_fiji_dist$Age_num_sd))#checking which rows contains NAs
## [1] 43

There in 1 NA in age withing data_fiji_dist.

data_fiji_dist<-data_fiji_dist[-c(43),]#eliminate row with NA
dim(data_fiji_dist)#number of observations. N=159 for Distributional equity
## [1] 159 141
str(data_fiji_dist$Distr_Eq_Impact)
##  chr [1:159] "4" "5" "5" "5" "4" "4" "4" "5" "4" "4" "3" "3" "5" "4" "4" ...
data_fiji_dist$Distr_Eq_Impact_int<-as.integer(data_fiji_dist$Distr_Eq_Impact)
str(data_fiji_dist$Distr_Eq_Impact_int)
##  int [1:159] 4 5 5 5 4 4 4 5 4 4 ...

PROCEDURAL EQUITY

ggplot(data_fiji, aes(Proced_Equity))+geom_bar()

Eliminate NA and don’t knows (including 6), and noanswer

data_fiji_proced <- data_fiji[data_fiji$Proced_Equity!=6,]

data_fiji_proced <- data_fiji_proced[data_fiji_proced$Proced_Equity!="dontknow",]

data_fiji_proced<- data_fiji_proced[data_fiji_proced$Proced_Equity!="noanswer",]

which(is.na(data_fiji_proced$Proced_Equity))
## integer(0)
data_fiji_proced<- data_fiji_proced[data_fiji_proced$Proced_Equity!="na",]
ggplot(data_fiji_proced, aes(Proced_Equity))+geom_bar()

dim(data_fiji_proced)#N=168 for Procedural Equity
## [1] 168 141
str(data_fiji_proced$Proced_Equity)
##  chr [1:168] "4" "5" "5" "5" "5" "4" "4" "4" "4" "3" "5" "4" "4" "4" "4" ...
data_fiji_proced$Proced_Equity_int<-as.integer(data_fiji_proced$Proced_Equity)
str(data_fiji_proced$Proced_Equity_int)
##  int [1:168] 4 5 5 5 5 4 4 4 4 3 ...

DISTRIBUTIONAL EQUITY

Reference categories: We will chose those with more responses (migrant1, females, married, secondary education). We only have to change the reference category for Education, because the other ones are already set.

data_fiji_dist<- within(data_fiji_dist, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts
data_fiji_dist$Distr_Eq_int<-data_fiji_dist$Distr_Eq_Impact_int
data_fiji_dist$Distr_Eq_int_cat<-as.factor(data_fiji_dist$Distr_Eq_int)

Bayesian (brms package)

The cumulative model assumes that the effect of a predictor is equal across the rating categories (e.g. the effect of marital status is the same across the levels of distributional equity). However, this assumption may not be appropriate, and marital status may impact rating categories differently. If the effects of predictors vary across rating categories, we call the resulting model to have category specific (cs) effects.

Checking proportional odds assumption with brms

(https://discourse.mc-stan.org/t/proportional-odds-assumptions/5316)

We need to investigate if each predictor has category specific effects. In order to do this we use an adjacent category model which uses the family=acat() instead of family=cumulative() (Burkner and Vuorre xx Ordinal regression models in psychology: a tutorial https://psyarxiv.com/x8swp/).

Additive model

#M.de5_add<-brm(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ #(1|village),data=data_fiji_dist,family=cumulative("logit"),control = list(adapt_delta = 0.99))  
#saveRDS(M.de5_add,file="M.de5_add.rda")
M.de5_add<-readRDS(file="M.de5_add.rda")
prior_summary(M.de5_add)# check default priors
##                 prior     class            coef   group resp dpar nlpar lb ub
##                (flat)         b                                              
##                (flat)         b      Age_num_sd                              
##                (flat)         b Education_cat31                              
##                (flat)         b Education_cat33                              
##                (flat)         b         Genderm                              
##                (flat)         b Marital_statuss                              
##                (flat)         b Marital_statusw                              
##                (flat)         b       migrant22                              
##                (flat)         b      wealth2_sd                              
##  student_t(3, 0, 2.5) Intercept                                              
##  student_t(3, 0, 2.5) Intercept               1                              
##  student_t(3, 0, 2.5) Intercept               2                              
##  student_t(3, 0, 2.5) Intercept               3                              
##  student_t(3, 0, 2.5) Intercept               4                              
##  student_t(3, 0, 2.5)        sd                                          0   
##  student_t(3, 0, 2.5)        sd                 village                  0   
##  student_t(3, 0, 2.5)        sd       Intercept village                  0   
##        source
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##  (vectorized)
##  (vectorized)
##  (vectorized)
##       default
##  (vectorized)
##  (vectorized)

Add weakly informative priors.

Weakly informative priors

priors2 <-prior(normal(0, 2.5), class = "b")+
  prior(normal(0, 2.5), class = "Intercept")+
  prior(student_t(3,0,2), class = "sd")

Modelling cs(Education_cat3) does not work.We will model it without specific categories for education

#```{r } M.de5_acat_all_noed<-brm(bf(Distr_Eq_int~cs(Gender)+cs(migrant2)+Education_cat3+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_acat_all_noed,file="M.de5_acat_all_noed.rda")
M.de5_acat_all_noed<-readRDS(file="M.de5_acat_all_noed.rda")
loo.M.de5_acat_all_noed<-loo(M.de5_acat_all_noed)
## Warning: Found 2 observations with a pareto_k > 0.7 in model
## 'M.de5_acat_all_noed'. It is recommended to set 'moment_match = TRUE' in order
## to perform moment matching for problematic observations.
M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")
loo.M.de5_acat_all2<-loo(M.de5_acat_all2)

#```{r } M.de5_acat_all2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_acat_all2,file="M.de5_acat_all2.rda")
M.de5_acat_all2<-readRDS(file="M.de5_acat_all2.rda")
loo.M.de5_acat_all2<-loo(M.de5_acat_all2)

Additive model

#```{r additive model2, cashe=TRUE} M.de5_add_2<-brm(bf(Distr_Eq_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_add_2,file="M.de5_add_2.rda") #

M.de5_add_2<-readRDS(file="M.de5_add_2.rda")
pars <- M.de5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_add_2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_add_2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_add_2%>%mcmc_trace

Chain convergence:

M.de5_add_2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_add_2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_add_2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_add_2<- as.array(M.de5_add_2)

Residual plot:

preds<-posterior_predict(M.de5_add_2,ndraws=250,summary=FALSE)

M.de5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_add_2.resids) #residuals look good

summary(M.de5_add_2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.26     0.02     1.02 1.00     2623     3148
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]       -4.31      0.58    -5.48    -3.24 1.00     3213     3217
## Intercept[2]       -2.40      0.44    -3.31    -1.55 1.00     3150     3174
## Intercept[3]       -1.82      0.43    -2.69    -1.02 1.00     3209     3027
## Intercept[4]        0.57      0.40    -0.22     1.37 1.00     3105     3089
## Genderm             0.09      0.37    -0.64     0.82 1.00     2960     3180
## migrant22          -0.71      0.35    -1.39    -0.03 1.00     3277     3008
## Education_cat31    -0.85      0.42    -1.66    -0.04 1.00     3208     3173
## Education_cat33    -0.99      0.40    -1.78    -0.19 1.00     3135     3093
## wealth2_sd          0.17      0.34    -0.52     0.82 1.00     3069     2976
## Age_num_sd          0.50      0.39    -0.26     1.29 1.00     3307     3174
## Marital_statuss     0.09      0.52    -0.97     1.11 1.00     3246     3172
## Marital_statusw    -1.33      0.52    -2.38    -0.32 1.00     3052     3071
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_add_2<-loo(M.de5_add_2)
loo.M.de5_add_2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -212.4  9.8
## p_loo        16.3  1.3
## looic       424.9 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare cumulative and acat

loo_compare(loo.M.de5_add_2,loo.M.de5_acat_all2,loo.M.de5_acat_all_noed)
##                     elpd_diff se_diff
## M.de5_add_2           0.0       0.0  
## M.de5_acat_all2      -1.3       1.5  
## M.de5_acat_all_noed -11.6       4.9

Gender x Migrant

#```{r genderxmigrant model2, cashe=TRUE} M.de5_int.mig2<-brm(bf(Distr_Eq_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_int.mig2,file="M.de5_int.mig2.rda") #

M.de5_int.mig2<-readRDS(file="M.de5_int.mig2.rda")
pars <- M.de5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.mig2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.mig2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.mig2%>%mcmc_trace

Chain convergence:

M.de5_int.mig2 %>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.mig2 %>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.mig2 %>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.mig2<- as.array(M.de5_int.mig2)
mcmc_acf(posterior_M.de5_int.mig2, 
         pars = c("b_Intercept[1]","b_Intercept[2]","b_Intercept[3]","b_Intercept[4]", "b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))

Residual plot:

preds<-posterior_predict(M.de5_int.mig2,ndraws=250,summary=FALSE)

M.de5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.mig2.resids) #residuals look good

summary(M.de5_int.mig2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.43      0.26     0.03     1.04 1.00     2558     2909
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -4.25      0.59    -5.43    -3.16 1.00     3084     3052
## Intercept[2]         -2.30      0.45    -3.22    -1.44 1.00     3215     3197
## Intercept[3]         -1.70      0.43    -2.58    -0.87 1.00     3205     2963
## Intercept[4]          0.73      0.42    -0.10     1.56 1.00     3339     3116
## Genderm               0.39      0.42    -0.45     1.22 1.00     3209     3131
## migrant22            -0.27      0.43    -1.13     0.59 1.00     2985     3161
## Education_cat31      -0.90      0.42    -1.71    -0.07 1.00     3102     3175
## Education_cat33      -1.06      0.39    -1.86    -0.32 1.00     2924     3183
## wealth2_sd            0.20      0.33    -0.45     0.86 1.00     3138     3091
## Age_num_sd            0.59      0.40    -0.19     1.37 1.00     3188     3120
## Marital_statuss       0.13      0.52    -0.88     1.16 1.00     3275     3153
## Marital_statusw      -1.34      0.54    -2.40    -0.29 1.00     3312     3127
## Genderm:migrant22    -1.19      0.70    -2.55     0.18 1.00     3112     3299
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
mcmc_intervals(posterior_M.de5_int.mig2, prob = 0.8, # 80% intervals
  prob_outer = 0.95, # 95%
  point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw","b_Genderm:migrant22"))

Conditional effects:

Ploting the interaction

conditions <- data.frame(Gender = c("f","m"))
conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects

inter_plot<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.95,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects
plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "darkred", 
                             "2"= "red3", 
                             "3"= "grey61",
                             "4" = "steelblue2",
                             "5" = "royalblue4")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))+
  scale_fill_discrete(labels = c("Very unfair", "Unfair", "Neutral","Fair", "Very fair"))
plot1

inter_plot.80<-conditional_effects(M.de5_int.mig2, "migrant2", prob=0.80,
                        categorical=TRUE, conditions=conditions) #re_formula=NULL takes into account random effects
plot2<- plot(inter_plot.80, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "darkred", 
                             "2"= "red3", 
                             "3"= "grey61",
                             "4" = "steelblue2",
                             "5" = "royalblue4")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))
plot2

Change colors

plot1

plot1<- plot(inter_plot, plot = FALSE)[[1]] + ##ggplot object
  theme_classic()+scale_color_manual(values = c("1"= "#CC0033", 
                             "2"= "#FF6699", 
                             "3"= "#CCCCCC",
                             "4" = "#CCCCFF",
                             "5" = "#6699FF")) +
  labs(x="Migrant Status")+
  theme(axis.text=element_text(size=10),
        axis.title=element_text(size=15),
        legend.text=element_text(size=10),
        legend.title=element_text(size=10))
plot1

loo.M.de5_int.mig2<-loo(M.de5_int.mig2)
loo.M.de5_int.mig2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -211.7 10.1
## p_loo        17.4  1.5
## looic       423.4 20.2
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x marital status

#{r M.de5_int.marital2, cashe=TRUE} M.de5_int.marital3<-brm(bf(Distr_Eq_int~Gender*Marital_status+Education_cat3+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior='yes', family=cumulative("logit"), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99)) #

#{r} saveRDS(M.de5_int.marital3,file="M.de5_int.marital3.rda") #

M.de5_int.marital3<-readRDS(file="M.de5_int.marital3.rda")
M.de5_int.marital2<-M.de5_int.marital3
pars <- M.de5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.marital2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.marital2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.marital2%>%mcmc_trace

Chain convergence:

M.de5_int.marital2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.marital2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.marital2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.marital2<- as.array(M.de5_int.marital2)

Residual plot:

preds<-posterior_predict(M.de5_int.marital2,ndraws=250,summary=FALSE)

M.de5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.marital2.resids) #residuals look good

summary(M.de5_int.marital2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.26     0.03     1.01 1.00     2371     2262
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -4.29      0.62    -5.57    -3.15 1.00     3287
## Intercept[2]               -2.37      0.47    -3.32    -1.47 1.00     3093
## Intercept[3]               -1.79      0.45    -2.70    -0.91 1.00     3180
## Intercept[4]                0.60      0.43    -0.21     1.46 1.00     3074
## Genderm                     0.10      0.42    -0.70     0.92 1.00     3175
## Marital_statuss             0.10      0.71    -1.30     1.47 1.00     2994
## Marital_statusw            -1.34      0.60    -2.53    -0.18 1.00     3119
## Education_cat31            -0.74      0.40    -1.54     0.03 1.00     3178
## Education_cat33            -0.97      0.39    -1.73    -0.20 1.00     3128
## wealth2_sd                  0.16      0.34    -0.51     0.83 1.00     2859
## Age_num_sd                  0.46      0.39    -0.29     1.23 1.00     3055
## migrant22                  -0.70      0.36    -1.39     0.00 1.00     3015
## Genderm:Marital_statuss    -0.05      0.95    -1.93     1.73 1.00     2974
## Genderm:Marital_statusw     0.08      0.92    -1.69     1.86 1.00     2526
##                         Tail_ESS
## Intercept[1]                2831
## Intercept[2]                3057
## Intercept[3]                3137
## Intercept[4]                2962
## Genderm                     3053
## Marital_statuss             3067
## Marital_statusw             2859
## Education_cat31             3023
## Education_cat33             3115
## wealth2_sd                  3027
## Age_num_sd                  3039
## migrant22                   2907
## Genderm:Marital_statuss     3158
## Genderm:Marital_statusw     3081
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.marital2<-loo(M.de5_int.marital2)
loo.M.de5_int.marital2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -215.0  9.8
## p_loo        18.4  1.6
## looic       430.1 19.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x education

#```{r M.de5_int.edu2, cashe=TRUE} M.de5_int.edu2<-brm(bf(Distr_Eq_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_int.edu2,file="M.de5_int.edu2.rda")
M.de5_int.edu2<-readRDS(file="M.de5_int.edu2.rda")
pars <- M.de5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.edu2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.edu2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.edu2%>%mcmc_trace

Chain convergence:

M.de5_int.edu2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.edu2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.edu2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.edu2<- as.array(M.de5_int.edu2)

Residual plot:

preds<-posterior_predict(M.de5_int.edu2,ndraws=250,summary=FALSE)

M.de5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.edu2.resids) #residuals look good

summary(M.de5_int.edu2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     0.98 1.00     2691     2685
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -4.28      0.59    -5.46    -3.17 1.00     3143
## Intercept[2]               -2.37      0.45    -3.26    -1.50 1.00     2702
## Intercept[3]               -1.79      0.43    -2.64    -0.95 1.00     2868
## Intercept[4]                0.60      0.42    -0.20     1.43 1.00     2926
## Genderm                     0.09      0.45    -0.77     0.99 1.00     3296
## Education_cat31            -0.84      0.58    -1.98     0.33 1.00     3286
## Education_cat33            -0.87      0.53    -1.89     0.13 1.00     3036
## Marital_statuss             0.09      0.54    -0.94     1.15 1.00     3022
## Marital_statusw            -1.30      0.54    -2.37    -0.26 1.00     3154
## wealth2_sd                  0.18      0.34    -0.51     0.87 1.00     3067
## Age_num_sd                  0.48      0.40    -0.31     1.26 1.00     3496
## migrant22                  -0.72      0.35    -1.40    -0.03 1.00     3202
## Genderm:Education_cat31     0.20      0.77    -1.31     1.72 1.00     3180
## Genderm:Education_cat33    -0.18      0.74    -1.58     1.24 1.00     3042
##                         Tail_ESS
## Intercept[1]                3303
## Intercept[2]                2992
## Intercept[3]                2988
## Intercept[4]                2848
## Genderm                     3217
## Education_cat31             3221
## Education_cat33             3258
## Marital_statuss             2997
## Marital_statusw             3194
## wealth2_sd                  2969
## Age_num_sd                  2895
## migrant22                   3174
## Genderm:Education_cat31     3112
## Genderm:Education_cat33     3010
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.edu2<-loo(M.de5_int.edu2)
loo.M.de5_int.edu2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -215.1  9.9
## p_loo        18.5  1.5
## looic       430.2 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x age

#```{r M.de5_int.age2, cashe=TRUE} M.de5_int.age2<-brm(bf(Distr_Eq_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.de5_int.age2,file="M.de5_int.age2.rda") #

M.de5_int.age2<-readRDS(file="M.de5_int.age2.rda")
pars <- M.de5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.age2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.age2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.age2%>%mcmc_trace

Chain convergence:

M.de5_int.age2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.age2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.age2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.age2<- as.array(M.de5_int.age2)

Residual plot:

preds<-posterior_predict(M.de5_int.age2,ndraws=250,summary=FALSE)

M.de5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.age2.resids) #residuals look good

summary(M.de5_int.age2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.26     0.03     1.03 1.00     2807     3058
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -4.31      0.58    -5.46    -3.22 1.00     2966     3173
## Intercept[2]          -2.40      0.44    -3.27    -1.54 1.00     3264     3144
## Intercept[3]          -1.81      0.43    -2.67    -0.99 1.00     3191     3061
## Intercept[4]           0.58      0.41    -0.23     1.40 1.00     3319     2529
## Genderm                0.09      0.37    -0.60     0.84 1.00     3289     3127
## Age_num_sd             0.45      0.57    -0.64     1.55 1.00     2917     2972
## Education_cat31       -0.85      0.42    -1.67    -0.03 1.00     2899     3215
## Education_cat33       -0.99      0.40    -1.78    -0.21 1.00     3189     3116
## Marital_statuss        0.08      0.52    -0.96     1.07 1.00     2873     3089
## Marital_statusw       -1.31      0.55    -2.39    -0.23 1.00     3052     2710
## wealth2_sd             0.17      0.35    -0.51     0.85 1.00     3329     2982
## migrant22             -0.72      0.36    -1.44    -0.04 1.00     3363     3172
## Genderm:Age_num_sd     0.08      0.66    -1.21     1.39 1.00     2693     3227
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.age2<-loo(M.de5_int.age2)
loo.M.de5_int.age2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -213.4  9.9
## p_loo        17.2  1.4
## looic       426.7 19.8
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x wealth

#```{r M.de5_int.wealth2, cashe=TRUE} M.de5_int.wealth2<-brm(bf(Distr_Eq_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_dist, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.de5_int.wealth2,file="M.de5_int.wealth2.rda")
M.de5_int.wealth2<-readRDS(file="M.de5_int.wealth2.rda")
pars <- M.de5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.de5_int.wealth2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.de5_int.wealth2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_dist$Distr_Eq_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.de5_int.wealth2%>%mcmc_trace

Chain convergence:

M.de5_int.wealth2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.de5_int.wealth2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.de5_int.wealth2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.de5_int.wealth2<- as.array(M.de5_int.wealth2)

Residual plot:

preds<-posterior_predict(M.de5_int.wealth2,ndraws=250,summary=FALSE)

M.de5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_dist$Distr_Eq_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.de5_int.wealth2.resids) #residuals look good

summary(M.de5_int.wealth2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Distr_Eq_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village) 
##    Data: data_fiji_dist (Number of observations: 159) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.37      0.25     0.02     0.97 1.00     2687     2897
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -4.35      0.59    -5.61    -3.24 1.00     3063     2965
## Intercept[2]          -2.42      0.44    -3.30    -1.58 1.00     2962     3015
## Intercept[3]          -1.84      0.42    -2.70    -1.02 1.00     2902     2824
## Intercept[4]           0.56      0.40    -0.21     1.33 1.00     2805     3181
## Genderm                0.08      0.36    -0.62     0.77 1.00     3057     3107
## wealth2_sd             0.41      0.41    -0.38     1.21 1.00     2749     2723
## Age_num_sd             0.51      0.39    -0.25     1.27 1.00     2864     3173
## Education_cat31       -0.84      0.41    -1.66    -0.06 1.00     3220     2848
## Education_cat33       -0.97      0.39    -1.76    -0.23 1.00     2974     3117
## Marital_statuss        0.06      0.53    -0.98     1.12 1.00     3056     3303
## Marital_statusw       -1.37      0.53    -2.39    -0.35 1.00     2755     2987
## migrant22             -0.71      0.35    -1.39    -0.04 1.00     3083     3134
## Genderm:wealth2_sd    -0.63      0.65    -1.91     0.68 1.00     3139     2822
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.de5_int.wealth2<-loo(M.de5_int.wealth2)
loo.M.de5_int.wealth2
## 
## Computed from 3200 by 159 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -212.6  9.9
## p_loo        16.6  1.4
## looic       425.3 19.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo.M.de5_int.wealth2,loo.M.de5_int.age2,loo.M.de5_int.edu2,loo.M.de5_int.marital2,loo.M.de5_int.mig2,loo.M.de5_add_2 )
##                    elpd_diff se_diff
## M.de5_int.mig2      0.0       0.0   
## M.de5_add_2        -0.8       2.0   
## M.de5_int.wealth2  -0.9       2.1   
## M.de5_int.age2     -1.7       2.0   
## M.de5_int.marital2 -3.3       2.1   
## M.de5_int.edu2     -3.4       2.1

PROCEDURAL EQUITY

data_fiji_proced<- within(data_fiji_proced, Education_cat3 <- relevel(Education_cat3, ref = 2)) #set secondary education as the reference category as it has more counts
str(data_fiji_proced$Proced_Equity_int)
##  int [1:168] 4 5 5 5 5 4 4 4 4 3 ...

#```{r } M.pe5_acat_all_noed<-brm(bf(Proced_Equity_int~cs(Gender)+cs(migrant2)+cs(Education_cat3)+cs(wealth2_sd)+cs(Age_num_sd)+cs(Marital_status)+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_acat_all_noed,file="M.pe5_acat_all_noed.rda")
M.pe5_acat_all_noed<-readRDS(file="M.pe5_acat_all_noed.rda")
M.pe5_acat_all_noed<-loo(M.pe5_acat_all_noed)

#```{r } M.pe5_acat_all2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=acat(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_acat_all2,file="M.pe5_acat_all2.rda")
M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")
M.pe5_acat_all2<-readRDS(file="M.pe5_acat_all2.rda")
loo.M.pe5_acat_all2<-loo(M.pe5_acat_all2)

Additive model

#```{r proced additive model2, cashe=TRUE} M.pe5_add_2<-brm(bf(Proced_Equity_int~Gender+migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_add_2,file="M.pe5_add_2.rda")
M.pe5_add_2<-readRDS(file="M.pe5_add_2.rda")
pars <- M.pe5_add_2%>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_add_2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_add_2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_add_2%>%mcmc_trace

Chain convergence:

M.pe5_add_2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_add_2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_add_2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_add_2<- as.array(M.pe5_add_2)

Residual plot:

preds<-posterior_predict(M.pe5_add_2,ndraws=250,summary=FALSE)

M.pe5_add_2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_add_2.resids) #residuals look good

summary(M.pe5_add_2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender + migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.42      0.27     0.02     1.02 1.00     2559     2732
## 
## Population-Level Effects: 
##                 Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]       -3.03      0.47    -3.99    -2.16 1.00     3324     3252
## Intercept[2]       -1.86      0.41    -2.70    -1.07 1.00     3356     3172
## Intercept[3]       -1.31      0.40    -2.11    -0.53 1.00     3304     3188
## Intercept[4]        0.63      0.39    -0.11     1.40 1.00     3351     3258
## Genderm             0.31      0.36    -0.40     1.00 1.00     3175     2941
## migrant22          -0.52      0.32    -1.14     0.11 1.00     3435     3256
## Education_cat31    -0.24      0.40    -1.03     0.55 1.00     3234     3161
## Education_cat33    -0.91      0.37    -1.62    -0.19 1.00     3234     3128
## wealth2_sd          0.10      0.31    -0.52     0.71 1.00     2907     2595
## Age_num_sd          0.21      0.35    -0.48     0.89 1.00     3331     3130
## Marital_statuss     0.55      0.50    -0.41     1.52 1.00     3228     3070
## Marital_statusw    -0.23      0.50    -1.15     0.76 1.00     3343     3064
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_add_2<-loo(M.pe5_add_2)
loo.M.pe5_add_2
## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -239.0  9.2
## p_loo        16.1  1.1
## looic       478.1 18.5
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare cumulative and acat

loo_compare(loo.M.pe5_add_2,loo.M.pe5_acat_all2,M.pe5_acat_all_noed)
##                     elpd_diff se_diff
## M.pe5_add_2           0.0       0.0  
## M.pe5_acat_all2      -0.9       1.6  
## M.pe5_acat_all_noed -11.8       5.6
posterior_M.pe5_add_2<- as.array(M.pe5_add_2)
mcmc_intervals(posterior_M.pe5_add_2, prob = 0.8, # 80% intervals
  prob_outer = 0.95, # 95%
  point_est = "mean"
,pars = c("b_Genderm", "b_migrant22", "b_Education_cat31","b_Education_cat33", "b_wealth2_sd","b_Age_num_sd","b_Marital_statuss","b_Marital_statusw"))

Gender x Migrant

#```{r proced genderxmigrant model2, cashe=TRUE} M.pe5_int.mig2<-brm(bf(Proced_Equity_int~Gender*migrant2+Education_cat3+wealth2_sd+Age_num_sd+Marital_status+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.mig2,file="M.pe5_int.mig2.rda")
M.pe5_int.mig2<-readRDS(file="M.pe5_int.mig2.rda")
pars <- M.pe5_int.mig2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.mig2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.mig2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.mig2%>%mcmc_trace

Chain convergence:

M.pe5_int.mig2 %>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.mig2 %>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.mig2 %>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

Residual plot:

preds<-posterior_predict(M.pe5_int.mig2,ndraws=250,summary=FALSE)

M.pe5_int.mig2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.mig2.resids) #residuals look good

summary(M.pe5_int.mig2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * migrant2 + Education_cat3 + wealth2_sd + Age_num_sd + Marital_status + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.27     0.02     1.03 1.00     2280     2620
## 
## Population-Level Effects: 
##                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]         -2.97      0.49    -3.96    -2.03 1.00     3280     3171
## Intercept[2]         -1.81      0.43    -2.68    -0.97 1.00     3294     3165
## Intercept[3]         -1.25      0.41    -2.06    -0.41 1.00     3206     3249
## Intercept[4]          0.69      0.41    -0.09     1.51 1.00     3259     3175
## Genderm               0.37      0.40    -0.41     1.17 1.00     3230     2941
## migrant22            -0.45      0.39    -1.22     0.30 1.00     3496     2510
## Education_cat31      -0.16      0.41    -0.96     0.62 1.00     3043     2908
## Education_cat33      -0.88      0.37    -1.63    -0.16 1.00     3163     3091
## wealth2_sd            0.12      0.32    -0.50     0.72 1.00     2914     3030
## Age_num_sd            0.20      0.37    -0.52     0.92 1.00     3306     2933
## Marital_statuss       0.57      0.50    -0.37     1.57 1.00     3259     3256
## Marital_statusw      -0.23      0.50    -1.17     0.75 1.00     3105     2972
## Genderm:migrant22    -0.17      0.67    -1.47     1.15 1.00     3210     3285
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.mig2<-loo(M.pe5_int.mig2)

Gender x marital status

#```{r M.pe5_int.marital2, cashe=TRUE} M.pe5_int.marital3<-brm(bf(Proced_Equity_int~Gender*Marital_status+Education_cat3+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#{r} saveRDS(M.pe5_int.marital3,file="M.pe5_int.marital3.rda") #

M.pe5_int.marital3<-readRDS(file="M.pe5_int.marital3.rda")
M.pe5_int.marital2<-M.pe5_int.marital3
pars <- M.pe5_int.marital2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.marital2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.marital2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <-  data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.marital2%>%mcmc_trace

Chain convergence:

M.pe5_int.marital2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.marital2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.marital2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

Residual plot:

preds<-posterior_predict(M.pe5_int.marital2,ndraws=250,summary=FALSE)

M.pe5_int.marital2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.marital2.resids) #residuals look good

summary(M.pe5_int.marital2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Marital_status + Education_cat3 + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     1.04 1.00     2905     2742
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -3.09      0.49    -4.10    -2.17 1.00     3157
## Intercept[2]               -1.92      0.43    -2.78    -1.10 1.00     3093
## Intercept[3]               -1.37      0.42    -2.19    -0.57 1.00     3109
## Intercept[4]                0.58      0.41    -0.20     1.38 1.00     3226
## Genderm                     0.14      0.40    -0.65     0.92 1.00     3327
## Marital_statuss             0.33      0.69    -1.00     1.68 1.00     3371
## Marital_statusw            -0.48      0.56    -1.55     0.62 1.00     3262
## Education_cat31            -0.11      0.41    -0.94     0.70 1.00     3240
## Education_cat33            -0.91      0.36    -1.64    -0.21 1.00     3110
## wealth2_sd                  0.15      0.31    -0.47     0.79 1.00     2986
## Age_num_sd                  0.22      0.36    -0.48     0.94 1.00     3421
## migrant22                  -0.54      0.34    -1.21     0.11 1.00     3176
## Genderm:Marital_statuss     0.50      0.92    -1.30     2.31 1.00     3429
## Genderm:Marital_statusw     1.08      1.04    -0.90     3.14 1.00     3298
##                         Tail_ESS
## Intercept[1]                3194
## Intercept[2]                3024
## Intercept[3]                3162
## Intercept[4]                3090
## Genderm                     3243
## Marital_statuss             3064
## Marital_statusw             3047
## Education_cat31             3169
## Education_cat33             3213
## wealth2_sd                  2806
## Age_num_sd                  3056
## migrant22                   3296
## Genderm:Marital_statuss     2980
## Genderm:Marital_statusw     3286
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.marital2<-loo(M.pe5_int.marital2)
loo.M.pe5_int.marital2
## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.5  9.2
## p_loo        17.8  1.2
## looic       480.9 18.4
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x education

#```{r M.pe5_int.edu2, cashe=TRUE} M.pe5_int.edu2<-brm(bf(Proced_Equity_int~Gender*Education_cat3+Marital_status+wealth2_sd+Age_num_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.edu2,file="M.pe5_int.edu2.rda")
M.pe5_int.edu2<-readRDS(file="M.pe5_int.edu2.rda")
pars <- M.pe5_int.edu2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.edu2%>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.edu2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.edu2%>%mcmc_trace

Chain convergence:

M.pe5_int.edu2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.edu2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.edu2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.edu2<- as.array(M.pe5_int.edu2)

Residual plot:

preds<-posterior_predict(M.pe5_int.edu2,ndraws=250,summary=FALSE)

M.pe5_int.edu2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.edu2.resids) #residuals look good

summary(M.pe5_int.edu2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Education_cat3 + Marital_status + wealth2_sd + Age_num_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.26     0.02     1.00 1.00     2776     2989
## 
## Population-Level Effects: 
##                         Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## Intercept[1]               -3.05      0.48    -4.01    -2.14 1.00     3201
## Intercept[2]               -1.89      0.43    -2.77    -1.05 1.00     3144
## Intercept[3]               -1.33      0.42    -2.16    -0.52 1.00     3236
## Intercept[4]                0.62      0.41    -0.16     1.43 1.00     3286
## Genderm                     0.24      0.44    -0.63     1.09 1.00     3284
## Education_cat31            -0.44      0.54    -1.51     0.60 1.00     3246
## Education_cat33            -0.90      0.45    -1.79    -0.02 1.00     3196
## Marital_statuss             0.53      0.50    -0.45     1.51 1.00     3352
## Marital_statusw            -0.16      0.51    -1.14     0.89 1.00     3253
## wealth2_sd                  0.08      0.32    -0.56     0.70 1.00     3169
## Age_num_sd                  0.21      0.36    -0.48     0.92 1.00     3069
## migrant22                  -0.53      0.33    -1.18     0.11 1.00     3303
## Genderm:Education_cat31     0.45      0.77    -1.06     1.96 1.00     3321
## Genderm:Education_cat33     0.01      0.70    -1.43     1.36 1.00     3329
##                         Tail_ESS
## Intercept[1]                2827
## Intercept[2]                3022
## Intercept[3]                3053
## Intercept[4]                3223
## Genderm                     3174
## Education_cat31             3031
## Education_cat33             3016
## Marital_statuss             3131
## Marital_statusw             3044
## wealth2_sd                  3171
## Age_num_sd                  3078
## migrant22                   3000
## Genderm:Education_cat31     3170
## Genderm:Education_cat33     3078
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.edu2<-loo(M.pe5_int.edu2)
loo.M.pe5_int.edu2
## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -241.1  9.3
## p_loo        18.4  1.2
## looic       482.2 18.7
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x age

#```{r M.pe5_int.age2, cashe=TRUE} M.pe5_int.age2<-brm(bf(Proced_Equity_int~Gender*Age_num_sd+Education_cat3+Marital_status+wealth2_sd+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#``

#saveRDS(M.pe5_int.age2,file="M.pe5_int.age2.rda")
M.pe5_int.age2<-readRDS(file="M.pe5_int.age2.rda")
pars <- M.pe5_int.age2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.age2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.age2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.age2%>%mcmc_trace

Chain convergence:

M.pe5_int.age2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.age2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.age2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.age2<- as.array(M.pe5_int.age2)

Residual plot:

preds<-posterior_predict(M.pe5_int.age2,ndraws=250,summary=FALSE)

M.pe5_int.age2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.age2.resids) #residuals look good

summary(M.pe5_int.age2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * Age_num_sd + Education_cat3 + Marital_status + wealth2_sd + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.41      0.27     0.03     1.05 1.00     2690     2872
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -3.04      0.48    -3.99    -2.12 1.00     3087     2954
## Intercept[2]          -1.87      0.42    -2.67    -1.05 1.00     2809     2971
## Intercept[3]          -1.32      0.40    -2.09    -0.51 1.00     2981     3131
## Intercept[4]           0.63      0.40    -0.12     1.38 1.00     2962     2877
## Genderm                0.30      0.36    -0.40     1.01 1.00     3196     3054
## Age_num_sd             0.26      0.50    -0.73     1.22 1.00     2974     3053
## Education_cat31       -0.24      0.41    -1.02     0.57 1.00     3342     3290
## Education_cat33       -0.91      0.36    -1.61    -0.21 1.00     3232     3017
## Marital_statuss        0.54      0.52    -0.45     1.59 1.00     3098     2827
## Marital_statusw       -0.25      0.53    -1.31     0.82 1.00     2922     2937
## wealth2_sd             0.11      0.31    -0.54     0.69 1.00     3280     3204
## migrant22             -0.52      0.32    -1.17     0.11 1.00     3290     3089
## Genderm:Age_num_sd    -0.08      0.64    -1.31     1.18 1.00     2900     2926
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.age2<-loo(M.pe5_int.age2)
loo.M.pe5_int.age2
## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.1  9.3
## p_loo        17.1  1.1
## looic       480.2 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Gender x wealth

#```{r M.de5_int.wealth2, cashe=TRUE} M.pe5_int.wealth2<-brm(bf(Proced_Equity_int~Gender*wealth2_sd+Age_num_sd+Education_cat3+Marital_status+migrant2+ (1|village)), data=data_fiji_proced, prior=priors2, sample_prior=‘yes’, family=cumulative(“logit”), iter=5000, warmup = 1000, chains = 4, cores=4, thin=5, refresh=0, control = list(adapt_delta = 0.99))

#```

#saveRDS(M.pe5_int.wealth2,file="M.pe5_int.wealth2.rda")
M.pe5_int.wealth2<-readRDS(file="M.pe5_int.wealth2.rda")
pars <- M.pe5_int.wealth2 %>% get_variables()
wch <- pars %>% stringr::str_detect("^b_((?!Intercept).)*$|^sd_.*") # Using a negative look-around
g <- purrr::map(pars[wch],
                ~(M.pe5_int.wealth2 %>%
                    hypothesis(paste(.x,"=0"), class="") %>%
                    plot(plot=FALSE))[[1]])
names(g) <- pars[wch]

patchwork::wrap_plots(g)

# g %>% sjPlot::plot_grid()

PP_check:

 yrep_char <- posterior_predict(M.pe5_int.wealth2)
 yrep_int <- sapply(data.frame(yrep_char, stringsAsFactors = TRUE), as.integer)
 y_int <- data_fiji_proced$Proced_Equity_int
ppc_bars(y_int, yrep_int)

Trace plots:

M.pe5_int.wealth2%>%mcmc_trace

Chain convergence:

M.pe5_int.wealth2%>% mcmc_plot(type='rhat_hist')
## Warning: Dropped 1 NAs from 'new_rhat(rhat)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

M.pe5_int.wealth2%>% mcmc_plot(type='neff_hist')
## Warning: Dropped 1 NAs from 'new_neff_ratio(ratio)'.
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Autocorrelation

M.pe5_int.wealth2%>% mcmc_plot(type='acf_bar')
## Warning: Removed 84 rows containing missing values (`geom_bar()`).

posterior_M.pe5_int.wealth2<- as.array(M.pe5_int.wealth2)

Residual plot:

preds<-posterior_predict(M.pe5_int.wealth2,ndraws=250,summary=FALSE)

M.pe5_int.wealth2.resids<-createDHARMa(simulatedResponse = t(preds),
                                 observedResponse = data_fiji_proced$Proced_Equity_int,
                                 fittedPredictedResponse = apply(preds,2,median),
                                 # 2 = we tell R to take preds and apply it tothe rows (1) or to the columns (2)
                                 integerResponse = FALSE)
plot(M.pe5_int.wealth2.resids) #residuals look good

summary(M.pe5_int.wealth2)
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: Proced_Equity_int ~ Gender * wealth2_sd + Age_num_sd + Education_cat3 + Marital_status + migrant2 + (1 | village) 
##    Data: data_fiji_proced (Number of observations: 168) 
##   Draws: 4 chains, each with iter = 5000; warmup = 1000; thin = 5;
##          total post-warmup draws = 3200
## 
## Group-Level Effects: 
## ~village (Number of levels: 10) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     0.40      0.27     0.02     1.03 1.00     2777     2744
## 
## Population-Level Effects: 
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]          -3.04      0.47    -3.97    -2.12 1.00     2865     2963
## Intercept[2]          -1.87      0.42    -2.70    -1.05 1.00     3069     3021
## Intercept[3]          -1.31      0.40    -2.10    -0.51 1.00     3137     2842
## Intercept[4]           0.63      0.40    -0.13     1.46 1.00     3284     3150
## Genderm                0.31      0.37    -0.41     1.04 1.00     3132     2800
## wealth2_sd             0.16      0.39    -0.61     0.94 1.00     2820     2527
## Age_num_sd             0.22      0.37    -0.51     0.95 1.00     3135     3091
## Education_cat31       -0.24      0.40    -1.03     0.57 1.00     3047     2839
## Education_cat33       -0.91      0.37    -1.65    -0.21 1.00     3074     2985
## Marital_statuss        0.55      0.52    -0.46     1.55 1.00     3234     2894
## Marital_statusw       -0.23      0.49    -1.21     0.70 1.00     3185     3064
## migrant22             -0.52      0.33    -1.16     0.12 1.00     3134     3241
## Genderm:wealth2_sd    -0.18      0.62    -1.40     1.03 1.00     3144     3004
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
loo.M.pe5_int.wealth2<-loo(M.pe5_int.wealth2)
loo.M.pe5_int.wealth2
## 
## Computed from 3200 by 168 log-likelihood matrix
## 
##          Estimate   SE
## elpd_loo   -240.1  9.3
## p_loo        17.0  1.1
## looic       480.1 18.6
## ------
## Monte Carlo SE of elpd_loo is 0.1.
## 
## All Pareto k estimates are good (k < 0.5).
## See help('pareto-k-diagnostic') for details.

Compare models

loo_compare(loo.M.pe5_int.wealth2,loo.M.pe5_int.age2,loo.M.pe5_int.edu2,loo.M.pe5_int.marital2,loo.M.pe5_int.mig2,loo.M.pe5_add_2 )
##                    elpd_diff se_diff
## M.pe5_add_2         0.0       0.0   
## M.pe5_int.wealth2  -1.0       0.3   
## M.pe5_int.age2     -1.1       0.2   
## M.pe5_int.mig2     -1.4       0.4   
## M.pe5_int.marital2 -1.4       1.1   
## M.pe5_int.edu2     -2.0       0.7